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Abstract 

Numerical simulations are a versatile tool providing insight into the complicated process of structure 
formation in cosmology [1]. This process is mainly governed by gravity, which is the dominant force 
on large scales. To date, a century after the formulation of general relativity [2], numerical codes 
for structure formation still employ Newton’s law of gravitation. This approximation relies on the 
two assumptions that gravitational fields are weak and that they are only sourced by non-relativistic 
matter. While the former appears well justified on cosmological scales, the latter imposes restrictions 
on the nature of the “dark” components of the Universe (dark matter and dark energy) which are, 
however, poorly understood. Here we present the first simulations of cosmic structure formation using 
equations consistently derived from general relativity. We study in detail the small relativistic effects 
for a standard ACDM cosmology which cannot be obtained within a purely Newtonian framework. Our 
particle-mesh N-body code computes all six degrees of freedom of the metric and consistently solves the 
geodesic equation for particles, taking into account the relativistic potentials and the frame-dragging 
force. This conceptually clean approach is very general and can be applied to various settings where 
the Newtonian approximation fails or becomes inaccurate, ranging from simulations of models with 
dynamical dark energy [3] or warm/hot dark matter [4] to core collapse supernova explosions [5]. 

The applicability of Newton’s law of gravitation in the context of cosmic structure formation has 
been discussed extensively in the recent literature [1 0 [5]. In particular, it is now well understood 
that this simplified description is quite accurate when applied within standard AGDM cosmology where 
perturbations come entirely from non-relativistic matter. However, the situation is not satisfactory for 
two reasons. Firstly, the quality of observational data is rapidly increasing, and upcoming galaxy surveys 
will eventually reach a precision where a naive treatment of the effects of spacetime geometry becomes 
insufficient PIIO]. Secondly, the true nature of dark matter and dark energy is not yet established. In order 
to study models beyond ACDM, some of which may feature relativistic sources of stress-energy (e.g. warm 
dark matter or dynamical dark energy), employing the Newtonian approximation is not always justified. 
A number of numerical codes has been developed for particular models [in na iia E], yet a general 
framework would be desirable. Furthermore, Newtonian gravity is acausal and not sensitive to the presence 
of a cosmological horizon. Even if a judicious interpretation of the output of Newtonian simulations can 
cure this problem at the linear level, it comes back when one goes beyond linear perturbation theory and 
it would be preferable to use the correct physics from the outset. 

Moving from the absolute space and time of the Newtonian picture towards a general relativistic view 
where geometry is dynamical poses a significant conceptual challenge. Recent progress is owed to a suitable 
formulation of the relativistic setting in terms of a weak-field expansion which is well adapted for (but not 
restricted to) cosmological applications O (TSl US]. Based on these ideas we have developed a numerical 
code, gevolution, which is designed to perform cosmological N-body simulations fully in the context of 
general relativity, evolving all six metric degrees of freedom. In brief, our approach can be summarized as 
follows. 

I. We choose a suitable ansatz for the metric which is split into background and perturbations. We 
work in Poisson gauge where the perturbed Friedmann-Lemaitre-Robertson-Walker metric is 

ds^ = a^ir) [— (1 -I- 24>) dr^ — 2Bidx^dT + (1 — 2$) Sijdx^dx^ + hijdx^dx^] , (I) 

where a denotes the scale factor of the background, a;* are comoving coordinates on the spacelike 
hypersurfaces, and r is conformal time. Bi (also denoted as B) is transverse and is responsible 
for frame dragging; hij is transverse and traceless, it contains the two spin-2 degrees of freedom of 
gravitational waves. 
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2. We assume that the perturbations of the metric, but not necessarily their derivatives, remain small on 
the scales of interest. This is a valid approximation for cosmological scales even when perturbations 
in the stress-energy tensor are large. 


3. Einstein’s equations are then expanded in terms of the metric perturbations but without a pertur¬ 
bative treatment of the stress-energy tensor. We include all terms linear in metric perturbations 
and go to quadratic order in terms which have two spatial derivatives acting on the scalar metric 
perturbations <i>, dr. This weak-field expansion contains all six components of the metric correctly 
at leading order. To determine the evolution of the metric, we numerically solve the “00” and the 
traceless part of the “ij” Einstein equations: 

(1 -t 4$) A$ - ^ (V$)^ = -'i-KGa^ST^ , (2a) 

-h”j + "H/iy - -^hij + j) + -!-($- dr) 

-b ($ - dr) . ($ - d')_^. -b 2 ($ - dr) ($ - d-) - 2 ($ - dr) d>,y -b 2$,,$,^ -b 

. (2b) 


= S^Ga" S,kTi - 
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Eq. (2a) is the generalisation of the Newtonian Poisson equation (Ad* = AnGa^Sp) and contains 
additional relativistic terms. Here Adr and 6Tq = —Sp = Tq — Tq are not required to be small - in 
fact, they become much larger than the background value Tq inside dense regions. Eq. (2b) can be 
used to evolve all the non-Newtonian degrees of freedom of the metric: d* — dr, the two spin-1 degrees 
of freedom, B, and the two spin-2 helicities, hij. To do this, we decompose the equation into scalar, 
vector and tensor parts in Fourier space. 


4. The metric is then used to solve the equations of motion for matter (and possibly other degrees of free¬ 
dom); collisionless particles propagate along geodesics in the perturbed geometry. This determines 
the evolution of the stress-energy tensor. 


The implementation in gevolution adopts a particle-mesh scheme where the metric field is represented 
on a regular cubic lattice while matter takes the form of an N-body ensemble of particles that samples 
the six-dimensional phase space. The stress-energy tensor on the lattice is obtained by a particle-to-mesh 
projection and is used in Einstein’s equations to solve for the metric perturbations. Finally, the particles 
are evolved by interpolating the metric to the particle positions and integrating the geodesic equation, 
gevolution is built on top of the LATfield2 c-b-b framework m- LATfield2 handles the lattice, the fields, 
the parallelization (using MPI) and the Fast Fourier Transforms which gevolution needs for solving Eqs. 
(2a) and (2b). We have extended the public version of LATfield2 to also handle the N-body particles, 
including projection and interpolation methods. We plan to make these features publicly available in the 
near future, together with a release of the gevolution code. 

As a first application we choose standard ACDM cosmology. We expect only small effects in this case, 
but we can use these simulations to gain confidence in our new approach. We generate [T5] two halo 
catalogs containing ^ 500,000 halos each, one using our relativistic approach and one with the traditional 
Newtonian approach as reference, starting from identical linear initial conditions. A Kolmogorov-Smirnov 
test shows no significant disagreement in the distributions of some 25 different halo properties such as 
mass, spin or shape parameters. 

The largest non-Newtonian effect is frame dragging, which is associated with the spin-1 perturbation 
B. This perturbation is sourced by the curl part of the momentum density found, e.g., in rotating massive 
objects. As long as perturbations are small this is a second-order effect which has been studied using 
perturbation theory [19] . The power spectrum for B has also been computed in the non-perturbative 
regime of structure formation using a post-Newtonian framework |20| which shares some features with 
our approach. These results are useful benchmarks for our code, but we now go beyond: our simulations 
track the full three-dimensional realisation of B (see Fig. [^. We can therefore measure the actual frame 
dragging force on individual particles. 

While the typical gravitational acceleration on Mpc scales is of the order of 10“® cm/s^, frame drag¬ 
ging contributes only at the level of 10“^^ cm/s^ (both numbers are mass-weighted rms values from our 
simulations at redshift zero), and the highest value we measure is some 10“^^ cm/s^. Thus, for objects 
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Figure 1: Spin-1 and spin-2 metric perturbations. Visualisation of a simulation volume of (512Mpc//i)® at 
redstiift z = 0. Dark matter halos are rendered as orange blobs scaled to the virial radius. Panel a illustrates the 
spin-2 perturbation hij by applying it as afhne transformation on spheroids of fixed radius; the shape and size of 
the resulting ellipsoids therefore indicates, respectively, the polarization and amplitude of the tensor perturbation. 
Panel b shows a stream plot of the spin-1 perturbation B, indicating how “spacetime is dragged around” by vortical 
matter flows. 


moving at 1000 km/s (a typical peculiar velocity at those scales) the Av due to frame dragging is no more 
than 1 km/s. It should be noted that these numbers are scale dependent, as the acceleration is larger on 
smaller scales. 

We compare our simulation results also with predictions [211122j from second order perturbation theory 
for the power spectra of hij, of B and of <i> — tk. The latter requires regularization in the infrared, which 
is implemented by the finite volume of the simulation. 

In order to fully capture the amplitude of the non-Newtonian terms, it is important that the scale of 
matter-radiation equality is represented in the simulation volume. On the other hand, these terms are 
generated by non-linearities and are therefore amplified at small scales. In order to obtain our results we 
had to cover at least three orders of magnitude in scale. Our largest simulation used a lattice with 4096^ 
points, corresponding to 6.7 x 10^° particles as we always start with one particle per grid point. 

Fig.j^shows the numerical power spectra at three redshifts. For these spectra we used eight simulations 
with lattices of 2048^ points and two different box sizes, all starting at redshift z = 100. Four simulations 
had a box size of (2048 Mpc//i)^, the other four were (512 Mpc//i)^. In the regime where it can be trusted 
we hnd excellent agreement with second order perturbation theory which demonstrates that our code 
produces valid results. At late times and on small scales, perturbation theory breaks down. Nonlinearities 
enhance the frame dragging by more than an order of magnitude at redshifts z = 1 and z = 0 and on 
scales k > 1/i/Mpc. The enhancement is even more dramatic for tensor perturbations. 

It should finally be noted that our framework jointly solves for background and perturbations in a 
self-consistent way. Therefore we confirm that clustering only has a small effect on the expansion rate of 
the Universe. The non-Newtonian effects which we have quantified in the nonlinear regime of structure 
formation remain small within a ACDM Universe, but they may nevertheless be measurable in the future. 
For the first time, general relativity has been implemented as the theory of gravity in a cosmological 
N-body code, making it possible to feed these effects back into the dynamics. Our numerical framework 
will be particularly useful in scenarios where relativistic sources are present, like in models of dynamical 
dark energy, topological defects, or with relativistic particles such as neutrinos or warm dark matter. Such 
scenarios are expected to display larger relativistic effects. Contrary to Newtonian schemes we are also able 
to solve the full geodesic equation for arbitrary velocities, allowing for a realistic propagation of radiation 
or high-velocity particles. 

Methods 

A weak-field expansion is useful in situations where the metric perturbation variables $, dt, Bi, hij defined 
in Eq. Q are small, as is the case on cosmological scales. In order to get a tractable set of equations, 
we retain all terms which are linear in metric perturbations, but from the quadratic terms we keep only 
the most relevant ones. These are built from the scalar perturbations $, 4^ only, and contain the highest 
possible number of spatial derivatives (two in this case, since all the partial differential equations are second 
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Figure 2: Power spectra. We show the power spectra of the gravitational potential $ (top, violet), the frame 
dragging potential B (second from top, blue), the difference of the relativistic potentials $ — (third from top, 
red) and the tensor perturbation hij (bottom, green) at redshifts z = 10, 1, 0. The black lines (dashed, dot- 
dashed, dot-dot-dashed, and dotted) indicate the corresponding results from second-order perturbation theory. As 
expected, the numerical results deviate from these perturbative extrapolations at low redshift and small scales. All 
perturbation variables are significantly enhanced on scales k > 1/i/Mpc at redshift z = 1 and below. Error bars 
indicate the random fluctuations from different realisations. 


order). This is enough to ensure that all six metric degrees of freedom are treated correctly at leading 
order, even in cases where some of them are strongly suppressed. 


Einstein’s equations. We determine the metric perturbations using the time-time as well as the trace¬ 
less space-space part of Einstein’s equations given in Eqs. (2a) and (2b). The remaining four equations 
are redundant, but we can use them as consistency check. The stress-energy tensor is constructed in the 
perturbed geometry and may hence also contain some terms which are linear in the metric perturbations 
- these terms have to be taken into account in order to maintain consistency. For a ACDM cosmology 
as the one studied here, the stress-energy tensor is obtained numerically by appropriate particle-to-mesh 
projections which are “dressed” by these geometric corrections. 

We then solve for the metric variables by treating the quadratic terms in Eqs. (2a) and (2b) as small 
perturbations. To this end, we simply estimate their value using the solutions taken from the previous 
time step and move them to the right-hand side. Since the equations are then approximated as linear in 
the metric perturbations, we can use Fourier methods to obtain the new solutions. This is convenient in 
particular since the gauge conditions can easily be implemented in Fourier space, where they reduce to 
local projection operations. One can check that the new solutions are stable by reinserting them into the 
quadratic terms and iterating the procedure. 


Particle trajectories. As dark matter is assumed to be collisionless, the particles move along geodesics. 
For non-relativistic velocities v = dx/dr the geodesic equation reads 

v'-fHv-bV^'-'HB-B'-bv X (V X B) =0, (3) 

where the B-dependent terms account for frame dragging. The integration of this equation can be simplified 
by writing v = v -|- B, which transforms Eq. ([^ to 

v' + nv + V^ + i^VBi = 0. (4) 

This new equation actually describes the acceleration as seen in a Gaussian ortho-normal coordinate 
frame. Frame dragging, characterized by the last term, has to compete with the “Newtonian” force V'l'. 
We measure both forces individually in our simulations. 
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Figure 3: Frame dragging. We show the frame dragging acceleration as a function of redshift z as measured 
in simulations of different spatial resolution. The solid lines correspond to mass-weighted rms values, whereas the 
dotted lines indicate the maximum value measured in the simulation volume, which was hxed to (512 Mpc//i)^. 
The acceleration is scale dependent and increases with resolution. Although not shown here, we checked that the 
effects of finite volume and realisation scatter are less relevant. 

Fig. i shows the mass-weighted rms value of the frame dragging acceleration (v • 9iB) for simulations 
with different spatial resolutions. The simulation volume was (512 Mpc/h)^ in all cases, but we used lattice 
sizes from 512^ to 4096^ points, reaching a best resolution of 125 kpc/h. The rms acceleration depends on 
the probed scale and therefore also on resolution. At our best resolution we remain slightly above galactic 
scales. At smaller scales we expect baryonic effects to become important. 

Power spectra. In Fig. we show the power spectra of the metric perturbations which are defined by 


47rfc^('I'(k, z)4'(k', z)) 

= (27r)35(3)(k-k')P4.(A;,z), 

(5a) 

Ank^{Bi(k, z)Bj{k', z)) 

= {2TT)^S^^'>{k-k/P,,PB{k,z), 

(5b) 

47rfc^(/iij(k, z)him{k',z)) 

= {2TT)^S^^Hk-k')M,pMk,z)- 

(5c) 


Here Pij = 6ij — kikj/k'^ is the transverse projector and the spin-2 projection operator is given by Mijim = 
PiiPjm + PimPji — PijPim- These dimensionless power spectra measure the amplitude square of the metric 
perturbations at scale k per log k interval. In order to suppress finite-volume and resolution effects we only 
measure scales which are at least 5 times smaller than the box size and have a wave number which is at 
least 5 times smaller than the Nyqvist frequency. 

Initial conditions. Initial data is generated using a linear input power spectrum at initial redshift 
(we start at z = 100) which can be obtained by running a Boltzmann code [23l [24] for the model. We 
use CLASS with the default cosmological parameters which describe a ACDM model (H^ = 0.0483, 
= He -l- Hf, = 0.312, our code currently has no special treatment for baryonic matter). 

Code availability. We will release the gevolution code, together with the necessary extensions of the 
LATfield2 library, on a public Git repository (https://github.eom/gevolution-code/gevolution-l.0.git). We 
strongly advise interested parties to wait for the fully documented version to become available. Prior to 
public release, access to the code can be given upon request, however, we do not offer extensive technical 
support in this case. 
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